##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(parallel)
library(ggplot2)
library(ggpubr)
library(cowplot)


##########################################################################################

option_list <- list(
  make_option(c("--OncoSG_rate_file_DGC"), type = "character") ,
  make_option(c("--OncoSG_rate_file_IGC"), type = "character") ,
  make_option(c("--njmu_rate_file"), type = "character") ,
  make_option(c("--tumor_type"), type = "character") ,
  make_option(c("--show_file"), type = "character") ,
  make_option(c("--out_path"), type = "character")
)

if(1!=1){
  
  tumor_type <- "STAD"
  OncoSG_rate_file_DGC <- "~/20220915_gastric_multiple/dna_combine/public_ref/OncoSG/OncoSG_STAD.mut_rate_DGC.tsv"
  OncoSG_rate_file_IGC <- "~/20220915_gastric_multiple/dna_combine/public_ref/OncoSG/OncoSG_STAD.mut_rate_IGC.tsv"
  njmu_rate_file <- "~/20220915_gastric_multiple/dna_combine/images/mutRate/MutRate.tsv"
  out_path <- ""
  show_file <- "~/20220915_gastric_multiple/dna_combine/mutsig_check/smg.list"
}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

OncoSG_rate_file_DGC <- opt$OncoSG_rate_file_DGC
OncoSG_rate_file_IGC <- opt$OncoSG_rate_file_IGC
njmu_rate_file <- opt$njmu_rate_file
tumor_type <- opt$tumor_type
show_file <- opt$show_file
out_path <- opt$out_path

###########################################################################################

oncosg_mut_rate_dgc <- fread(OncoSG_rate_file_DGC,sep = "\t",header = T)
oncosg_mut_rate_igc <- fread(OncoSG_rate_file_IGC,sep = "\t",header = T)
njmu_mut_rate <- fread(njmu_rate_file,sep = "\t",header = T)
show_gene <- data.frame(fread(show_file , header = F))$V1

###########################################################################################
## 提取SMG
oncosg_mut_rate_dgc$SampleNum <- 33
oncosg_mut_rate_dgc$Class <- "DGC"
oncosg_mut_rate_igc$SampleNum <- 63
oncosg_mut_rate_igc$Class <- "IGC"
oncosg_mut_rate <- rbind(oncosg_mut_rate_dgc,oncosg_mut_rate_igc)
njmu_mut_rate <- subset( njmu_mut_rate , Hugo_Symbol %in% show_gene )
oncosg_mut_rate <- subset( oncosg_mut_rate , gene %in% show_gene )

njmu_mut_rate <- njmu_mut_rate[,c(2,1,3:5)]
colnames(njmu_mut_rate)[1] <- "gene"


oncosg_mut_rate <- oncosg_mut_rate[,-c(4:7)]
colnames(oncosg_mut_rate) <- c("gene","MutNum","MutRate","SampleNum","Class")
oncosg_mut_rate <- oncosg_mut_rate[,c("gene","Class","MutNum","SampleNum","MutRate")]

#更改njmu_longtowide格式
library(reshape2)
class <- c("DGC","IGC")
njmu_mut_rate <- njmu_mut_rate[njmu_mut_rate$Class %in% class]
n <- colnames(njmu_mut_rate[,-c(1:2)])
longdata2 <- melt(njmu_mut_rate,id.vars=c("gene","Class"),measure.vars=n)
njmu_mut_rate <- dcast(longdata2,gene~variable+Class)

#njmu_mut_rate_dgc <- njmu_mut_rate[njmu_mut_rate$Class=="DGC"][,c(1,3:5)]
#njmu_mut_rate_igc <- njmu_mut_rate[njmu_mut_rate$Class=="IGC"][,c(1,3:5)]
#colnames(njmu_mut_rate_dgc)[2:ncol(njmu_mut_rate_dgc)] <- paste0( colnames(njmu_mut_rate_dgc)[2:ncol(njmu_mut_rate_dgc)] , "_DGC" )
#colnames(njmu_mut_rate_igc)[2:ncol(njmu_mut_rate_igc)] <- paste0( colnames(njmu_mut_rate_igc)[2:ncol(njmu_mut_rate_igc)] , "_IGC" )
#njmu_mut_rate <- merge(njmu_mut_rate_dgc, njmu_mut_rate_igc, by= "gene" ,all.y=T)

#更改oncosg_longtowide格式
n <- colnames(oncosg_mut_rate[,-c(1:2)])
longdata2 <- melt(oncosg_mut_rate,id.vars=c("gene","Class"),measure.vars=n)
oncosg_mut_rate <- dcast(longdata2,gene~variable+Class)


#更改后缀名称
colnames(njmu_mut_rate)[2:ncol(njmu_mut_rate)] <- paste0( colnames(njmu_mut_rate)[2:ncol(njmu_mut_rate)] , "_njmu" )
colnames(oncosg_mut_rate)[2:ncol(oncosg_mut_rate)] <- paste0( colnames(oncosg_mut_rate)[2:ncol(oncosg_mut_rate)] , "_oncosg" )


## 比较突变率
dat_rate <- merge(oncosg_mut_rate, njmu_mut_rate, by= "gene",all.y=T)


###########################################################################################


###########################################################################################
## 计算显著性
result1 <- Reduce(function(x,y)bind_rows(x,y),mclapply(dat_rate$gene,function(geneN){
  print(geneN)
  
  tmp <- subset(dat_rate , gene == geneN )
  
  ## DGC
  a <- tmp$MutNum_DGC_njmu
  c <- round(tmp$MutNum_DGC_njmu/tmp$MutRate_DGC_njmu) - tmp$MutNum_DGC_njmu
  b <- tmp$MutNum_DGC_oncosg
  d <- round(tmp$MutNum_DGC_oncosg/tmp$MutRate_DGC_oncosg) - tmp$MutNum_DGC_oncosg
  
  if(is.na(a)){a=0}
  if(is.na(b)){b=0}
  if(is.na(c)){c=0}
  if(is.na(d)){d=0}
  
  result=fisher.test(matrix(c(a,b,c,d),nrow=2))
  
  p=result[["p.value"]]
  OR=round(result[["estimate"]][["odds ratio"]],3)
  
  tmp$p_DGC <- p
  tmp$OR_DGC <- OR
  
  ## IGC
  a <- tmp$MutNum_IGC_njmu
  c <- round(tmp$MutNum_IGC_njmu/tmp$MutRate_IGC_njmu) - tmp$MutNum_IGC_njmu
  b <- tmp$MutNum_IGC_oncosg
  d <- round(tmp$MutNum_IGC_oncosg/tmp$MutRate_IGC_oncosg) - tmp$MutNum_IGC_oncosg
  
  if(is.na(a)){a=0}
  if(is.na(b)){b=0}
  if(is.na(c)){c=0}
  if(is.na(d)){d=0}
  
  result=fisher.test(matrix(c(a,b,c,d),nrow=2))
  
  p=result[["p.value"]]
  OR=round(result[["estimate"]][["odds ratio"]],3)
  
  tmp$p_IGC <- p
  tmp$OR_IGC <- OR
  
  tmp
  
},mc.cores=1))

result1 <- result1[order(result1$MutRate_IGC_njmu , decreasing = T),]
result1[is.na(result1)] <- 0


out_name <- paste0(out_path , "/OncoSG_" , tumor_type , ".mutRate.compare_select.tsv")
write.table(result1,out_name,sep="\t" , quote =F , row.names = F )

###########################################################################################

## 展示所有，上面为IGC，下面为DGC

result2 <- result1[,c("gene" , "MutRate_IGC_njmu" , "MutRate_IGC_oncosg" , "OR_IGC" , "p_IGC")]
result3 <- result1[,c("gene" , "MutRate_DGC_njmu" , "MutRate_DGC_oncosg" , "OR_DGC" , "p_DGC")]

result2$Type <- "IGC"
result3$Type <- "DGC"

colnames(result2) <- c("gene" , "njmu_rate" , "oncosg_rate" , "OR" , "p" , "Type")
colnames(result3) <- c("gene" , "njmu_rate" , "oncosg_rate" , "OR" , "p" , "Type")

result4 <- rbind( result2 , result3 )

###########################################################################################
## 柱状图展示
## p值
result_use <- result4
result_use$p_text=ifelse(result_use$p>=0.05,"","*")
result_use$p_text=ifelse(result_use$p<0.05 & result_use$p>0.01,"*",result_use$p_text)
result_use$p_text=ifelse(result_use$p<0.01 & result_use$p>0.001,"**",result_use$p_text)
result_use$p_text=ifelse(result_use$p<0.001 ,"***",result_use$p_text)

gene_order <- unique(result_use$gene)

##########################################################################################

result2 <- melt(result_use[,c("gene","njmu_rate","oncosg_rate","p_text","Type")])
result2$p_pos <- 0.9
result2$percent_pos <- result2$value + 0.003

result2$value_percent=paste(round(result2$value * 100),"%",sep="")
result2$gene <- factor( result2$gene , levels = gene_order , order = T )
result2$variable <- ifelse( result2$variable == "njmu_rate" , "NJMU" , "oncosg" )

result2$value[result2$Type == "DGC"] <- -result2$value[result2$Type == "DGC"]
result2$p_pos[result2$Type == "DGC"] <- -result2$p_pos[result2$Type == "DGC"]
result2$percent_pos[result2$Type == "DGC"] <- -result2$percent_pos[result2$Type == "DGC"] -0.04
result2 <- within(result2, Type <- factor(Type, levels = c("IGC", "DGC")))

p <- ggplot(result2,mapping = aes(x = gene , y = value , fill = variable)) + geom_bar(stat = 'identity', position = 'dodge') +
  facet_grid(vars(Type),scales = "free")+
  theme_bw() +
  theme(axis.text.y = element_blank())+
  theme(panel.grid=element_blank())+labs(x = 'Genes',y = 'Muatetation R') +
  theme(axis.title =element_text(size = 15),axis.text =element_text(size = 12, color = 'black'))+
  geom_text(aes(label= value_percent , x = gene, y= percent_pos ), position=position_dodge(0.9), vjust=0 ,color="black", face='bold' , family="Helvetica")+
  geom_text(aes(label=p_text , x = gene , y = p_pos ),size=7,family="Helvetica")+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1 , size = 15 , family="Helvetica"))+
  xlab("") +
  theme(
    title =element_text(size=4, face='bold'),
    legend.title = element_blank(),
    legend.text = element_text(size = 12),
    legend.key.width = unit(1, "cm"),
    legend.key.height = unit(1, "cm"),
    plot.title = element_text(size = 30, face = "bold"),
    text = element_text(family="Helvetica")
  )


width <- 10
height <- 7
images_name <- paste0(out_path , "/MutRate.compare_OncoSG.DriverGene.pdf")
ggsave( images_name , p , width = width , height = height )


